import pandas as pd
import numpy as np
import math
from scipy.special import factorial
from sympy import symbols, sympify, exp, N as high_precision_evaluator


epsilon = 0.15
N = 16
D = 4
R = N / D
sigma = 0.01
d = 1
lambda_ = 0.00015
Delta_d = lambda_ / 2
M = 40
K = 40
W = 1

# 假设卫星的坐标 (x_sat, y_sat, z_sat)
# 地球同步轨道卫星的大致高度为 35786 km，假设在 z 轴上，x 和 y 为 0
satellite_coords = (0, 0, 3578)


file1_path = "./cluster_info.xlsx"
file2_path = "special_square_points.xlsx"

# 获取分组数据
df1 = pd.read_excel(file1_path, sheet_name="Sheet1")
df2 = pd.read_excel(file2_path, sheet_name="Sheet1")

# 保存数据
grouped_data = {
    "Cluster_Label": df1["Cluster Label"],
    "Tar_Point_X": df1["Tar Point X"],
    "Tar_Point_Y": df1["Tar Point Y"],
    "ClusterCenterX": df1["ClusterCenterX"],
    "ClusterCenterY": df1["ClusterCenterY"],
    "Number_of_Covert_Points": df1["Number of Covert Points"],
    "DistanceToCovertClusterCenter": df1["DistanceToCovertClusterCenter"],
}
# 将group_data变为数组，每一项当中的元素为一个字典, 字典的键为"Cluster_Label", "Tar_Point_X", "Tar_Point_Y", "ClusterCenterX", "ClusterCenterY", "Number_of_Covert_Points", "DistanceToCovertClusterCenter"
group_data = []
for i in range(len(grouped_data["Cluster_Label"])):
    group_data.append(
        {
            "Cluster_Label": grouped_data["Cluster_Label"][i],
            "Tar_Point_X": grouped_data["Tar_Point_X"][i],
            "Tar_Point_Y": grouped_data["Tar_Point_Y"][i],
            "ClusterCenterX": grouped_data["ClusterCenterX"][i],
            "ClusterCenterY": grouped_data["ClusterCenterY"][i],
            "Number_of_Covert_Points": grouped_data["Number_of_Covert_Points"][i],
            "DistanceToCovertClusterCenter": grouped_data[
                "DistanceToCovertClusterCenter"
            ][i],
        }
    )

# 将group_data按照"Number_of_Covert_Points"排序
group_data.sort(key=lambda x: x["Number_of_Covert_Points"], reverse=True)
# print(group_data[6])
# print(group_data[3])
special_points = {
    "Special_Point_X": df2["X"],
    "Special_Point_Y": df2["Y"],
}
# 将special_points变为数组，每一项当中的元素为一个字典, 字典的键为"Special_Point_X", "Special_Point_Y"
special_points = []
for i in range(len(df2["X"])):
    special_points.append(
        {
            "Special_Point_X": df2["X"][i],
            "Special_Point_Y": df2["Y"][i],
        }
    )


#########
# 计算逻辑
def calculate_angles(satellite_coords, ground_node_coords):
    # 计算从卫星到地面节点的向量
    vector = np.array(satellite_coords) - np.array(ground_node_coords)
    # 计算向量的模长
    r = np.linalg.norm(vector)
    # 计算 theta 角度（与 z 轴的夹角）
    theta_ai = np.arccos(vector[2] / r)

    # 计算 phi 角度（在 x-y 平面的角度）
    phi_ai = np.arctan2(vector[1], vector[0])
    # 计算 theta 和 phi 的正弦值
    # sin_theta_ai = np.sin(theta_ai)
    # sin_phi_ai = np.sin(phi_ai)

    # 需要保证theta_ai和phi_ai是锐角
    if theta_ai > np.pi / 2:
        theta_ai = np.pi - theta_ai

    if phi_ai > np.pi / 2:
        phi_ai = np.pi - phi_ai

    return theta_ai, phi_ai


def calculate_hw(M, K, theta_b, theta_w, phi_b, phi_w, lambda_, Delta_d):
    # 计算垂直相位差 Theta
    Theta = (math.sin(theta_b) - math.sin(theta_w)) * math.pi * Delta_d / lambda_
    # 计算水平相位差 Phi
    Phi = (math.sin(phi_b) - math.sin(phi_w)) * math.pi * Delta_d / lambda_

    # 计算分子部分 sin(M*Theta) / sin(Theta) 的平方
    if math.isclose(math.sin(Theta), 0, abs_tol=1e-15):
        numerator_theta = M**2
    else:
        numerator_theta = (math.sin(M * Theta) / math.sin(Theta)) ** 2

    # 计算分子部分 sin(K*Phi) / sin(Phi) 的平方
    if math.isclose(math.sin(Phi), 0, abs_tol=1e-15):
        numerator_phi = K**2
    else:
        numerator_phi = (math.sin(K * Phi) / math.sin(Phi)) ** 2

    # 计算 |hw|² 的结果
    hw_squared = (1 / (M * K)) * numerator_theta * numerator_phi
    return hw_squared


def calculate_P_a(epsilon, N, sigma, d, alpha, hw):
    # 计算 P_a 的值
    P_a = (2 * epsilon * math.sqrt(N) * sigma) / (d ** (-alpha) * hw * N)
    return P_a


# 从group_data从提取出"Tar_Point_X", "Tar_Point_Y", "ClusterCenterX", "ClusterCenterY"形成坐标点，形成两个元组数组，第一个数组为Tar_Point_X, "Tar_Point_Y", 第二个数组为"ClusterCenterX", "ClusterCenterY"，所有的z坐标为0
# 计算角度
angles_tar = []
for i in range(len(group_data)):
    ground_node_coords = (
        group_data[i]["Tar_Point_X"],
        group_data[i]["Tar_Point_Y"],
        0,
    )
    theta_ai, phi_ai = calculate_angles(satellite_coords, ground_node_coords)
    angles_tar.append((theta_ai, phi_ai))

angles_center = []
for i in range(len(group_data)):
    ground_node_coords = (
        group_data[i]["ClusterCenterX"],
        group_data[i]["ClusterCenterY"],
        0,
    )
    theta_ai, phi_ai = calculate_angles(satellite_coords, ground_node_coords)
    angles_center.append((theta_ai, phi_ai))
# 遍历所有的angles_center如果发现里面的theta_ai或phi_ai的绝对值如果小于等于0.1，就将其变为0.3
for i in range(len(angles_center)):
    # if abs(angles_center[i][0]) <= 0.1:
    #     angles_center[i] = (0.3, angles_center[i][1])
    if abs(angles_center[i][1]) <= 0.1:
        angles_center[i] = (angles_center[i][0], 0.3)


# 输出angles_center[4]与angles_center[3]
# print("angles_center[6]", angles_center[6])
# print("angles_center[3]", angles_center[3])
# 处理special_points
angles_special = []
for i in range(len(special_points)):
    ground_node_coords = (
        special_points[i]["Special_Point_X"],
        special_points[i]["Special_Point_Y"],
        0,
    )
    theta_ai, phi_ai = calculate_angles(satellite_coords, ground_node_coords)
    angles_special.append((theta_ai, phi_ai))


# 根据angles_tar与angles_center计算hw
hw_tar = []
for i in range(len(angles_tar)):
    theta_b, phi_b = angles_tar[i]
    theta_w, phi_w = angles_center[i]
    hw = calculate_hw(M, K, theta_b, theta_w, phi_b, phi_w, lambda_, Delta_d)
    hw_tar.append(hw)

# 根据(0, 0, 0)与angles_center计算hw
hw_willie = []
for i in range(len(angles_center)):
    theta_b, phi_b = angles_center[i]
    theta_w, phi_w = (0, 0)
    hw = calculate_hw(M, K, theta_b, theta_w, phi_b, phi_w, lambda_, Delta_d)
    hw_willie.append(hw)
# 根据angles_special与angles_special计算hw
hw_special = []
for i in range(len(angles_special)):
    theta_b, phi_b = angles_special[i]
    theta_w, phi_w = angles_special[i]
    hw = calculate_hw(M, K, theta_b, theta_w, phi_b, phi_w, lambda_, Delta_d)
    hw_special.append(hw)

# 根据(0, 0, 0)与angles_special计算hw

hw_special_willie = []
for i in range(len(angles_special)):
    theta_b, phi_b = angles_special[i]
    theta_w, phi_w = (0, 0)
    hw = calculate_hw(M, K, theta_b, theta_w, phi_b, phi_w, lambda_, Delta_d)
    hw_special_willie.append(hw)

# 输出结果
print("hw_tar:", hw_tar)
print("hw_willie:", hw_willie)
print("hw_special:", hw_special)
print("hw_special_willie:", hw_special_willie)


# 计算满足willie处隐蔽约束的P_a，然后就能推导出真正的P_a，这段代码就是在计算真正的P_a
P_a_willie = []
for i in range(len(hw_willie)):
    P_a = calculate_P_a(epsilon, N, sigma, d, 0, hw_willie[i])
    P_a_willie.append(P_a)

print("P_a_willie:", P_a_willie)

# 根据hw_special_willie计算P_a
P_a_willie_special = []
for i in range(len(hw_special_willie)):
    P_a = calculate_P_a(epsilon, N, sigma, d, 0, hw_special_willie[i])
    P_a_willie_special.append(P_a)


# willie处的pa乘hw_tar为P_a_tar
P_a_tar = []
for i in range(len(P_a_willie)):
    P_a = P_a_willie[i] * hw_tar[i]
    P_a_tar.append(P_a)

# 根据P_a_willie_special与hw_special计算P_a_tar_special
P_a_tar_special = []
for i in range(len(P_a_willie_special)):
    P_a = P_a_willie_special[i] * hw_special[i]
    P_a_tar_special.append(P_a)


# 计算 beta 的值
beta = -np.sqrt(N / (2 * np.pi * (np.exp(2 * R) - 1)))

# 计算v的值
v = exp(R) - 1

# 注意下面的函数中的K为K_ri
K_ri = 3

# 这里P_b写死为0
P_b = 0
sigma_b = sigma
Omega = 0


def P_beCal(beta, v, K, sigma_b, Omega, P_b, P_a):
    x = symbols("x")
    # 使用有限的上界代替无穷级数
    upperBound = 5

    # 构建内部的双重求和表达式
    def inner_sum():
        total = 0
        for q in range(0, upperBound + 1):
            for p in range(0, upperBound + 1):
                term = (
                    (K**q)
                    * ((K * (1 + K) * (sigma_b**2 + Omega * P_b) * x / P_a) ** p)
                    / (factorial(p) * factorial(p + q))
                )
                total += term
        return total

    innerSum = high_precision_evaluator(inner_sum(), 64)
    # 被积函数
    integrand = (
        1 - exp(-(K + ((1 + K) * (sigma_b**2 + Omega * P_b) * x) / P_a)) * innerSum
    )
    # 将积分写成矩形近似的形式，区间划分为 100 份
    dx = (v - 1 / (2 * beta) - (v + 1 / (2 * beta))) / 100
    x_values = np.arange(v + 1 / (2 * beta), v - 1 / (2 * beta), dx)
    # 计算 innerExp
    # innerExp = np.exp(
    #     -(np.array([K]) + ((1 + K) * (sigma_b**2 + Omega * P_b) * x_values) / P_a)
    # )
    # 显示 innerExp 的信息
    # print(f"this is innerExp {innerExp}")
    # 计算积分近似
    integrand_values = np.array(
        [sympify(integrand.subs(x, val)).evalf() for val in x_values]
    )

    # Calculate the integral approximation
    integral_approx = dx * np.sum(integrand_values) * -beta
    return integral_approx


def calculate_min_AoI(N, W, P_be):
    """
    此函数用于计算 minAoI 的值
    :param N: 输入参数 N
    :param W: 输入参数 W
    :param P_be: 输入参数 P_be
    :return: 计算得到的 minAoI 的值
    """
    min_AoI = (2 * N) / (W * (1 - P_be)) + (N / (2 * W))
    return min_AoI


print("P_a_tar:", P_a_tar)

# 计算P_be
P_be = []
for i in range(len(P_a_tar)):
    P_be.append(P_beCal(beta, v, K_ri, sigma_b, Omega, P_b, P_a_tar[i]))
# 将P_a_tar_special的值也带入计算，加入到P_be中
for i in range(len(P_a_tar_special)):
    P_be.append(P_beCal(beta, v, K_ri, sigma_b, Omega, P_b, P_a_tar_special[i]))

# 计算min_AoI
min_AoI = []
for i in range(len(P_be)):
    min_AoI.append(calculate_min_AoI(N, W, P_be[i]))

print("P_be:", P_be)
# 将min_AoI的每一个值加上前面的min_AoI
for i in range(1, len(min_AoI)):
    min_AoI[i] = min_AoI[i] + min_AoI[i - 1]

AoI_Sum = 0
# 处理group_data
for i in range(len(group_data)):
    AoI_Sum += min_AoI[i] * group_data[i]["Number_of_Covert_Points"]
# 处理special_points，这里的索引应该从len(group_data)开始
for i in range(len(special_points)):
    AoI_Sum += min_AoI[i + len(group_data)]


# 对group_data中的"Number_of_Covert_Points"进行求和
Number_of_Covert_Points_Sum = 0
for i in range(len(group_data)):
    Number_of_Covert_Points_Sum = (
        Number_of_Covert_Points_Sum + group_data[i]["Number_of_Covert_Points"]
    )
# 对special_points中的"Number_of_Covert_Points"进行求和
for i in range(len(special_points)):
    Number_of_Covert_Points_Sum = Number_of_Covert_Points_Sum + 1


res = AoI_Sum / Number_of_Covert_Points_Sum

print("AoI_Sum:", AoI_Sum)
print("Number_of_Covert_Points_Sum:", Number_of_Covert_Points_Sum)
print(res)

# 如过res小于500，写入文件


if res < 900:
    # 再以追加模式打开文件添加数据
    with open("res.txt", "a") as f:
        f.write(str(res) + "\n")
